Introduction to R for RNA-Seq

Acknowledgements

This material is adapted from an R for RNA-Seq workshop originally run here.

R for RNA-Seq course

This course has been designed to introduce biologists to R for RNA-Seq analysis. The focus here is on using tidyverse to analyse RNA-Seq data, as we believe this is a productive and engaging way to learn R for RNA-Seq analysis. In this course we will use some new packages, ttBulk and tidyHeatmap. These packages provide a friendly tidyverse-style way to perform analysis of RNA-Seq data.

Setup

Data files

Data files are available from the data folder in GitHub here. You should download the files listed below and place them into a folder called data in your working directory.

Data files:

• GSE60450_GeneLevel_Raw_data.csv
• GSE60450_filtered_metadata.csv

R Packages

Packages used:

• tidyverse
• ttBulk
• tidyHeatmap
• edgeR
• devtools

Here we will use packages from the 3 main repositories of R packages: Bioconductor, CRAN and GitHub. To install the packages you can follow the steps below.

R version

• Install R >= 3.6. You must have R >= 3.6 to use ttBulk and tidyHeatmap.

Bioconductor packages

• Get the latest version of Bioconductor and edgeR package by starting R and entering the commands:

if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("edgeR")

CRAN packages

Install the CRAN packages with the command below.

install.packages(c("tidyverse", "devtools"))

GitHub packages

ttBulk will be added to Bioconductor and when it is added there you can install it using the usual Bioconductor commands. tidyHeatmap will be added to CRAN and when it is there you can install it with the install.packages() command. In the meantime you can install ttBulk and tidyHeatmap from their development sites in Github with the commands below.

# install ttBulk
devtools::install_github("stemangiola/ttBulk@dev")

# install tidyHeatmap
devtools::install_github("stemangiola/tidyHeatmap")

Overview

TODO: Maybe add workflow image showing steps and functions used

• Reading in tables of counts and metadata
• Filtering and Normalisation of counts
• Differential expression analysis
• Visualisation (Heatmaps, MA plot, Volcano plot)

Introduction and data import

Measuring gene expression on a genome-wide scale has become common practice over the last two decades or so, with microarrays predominantly used pre-2008. With the advent of next generation sequencing technology in 2008, an increasing number of scientists use this technology to measure and understand changes in gene expression in often complex systems. As sequencing costs have decreased, using RNA-Seq to simultaneously measure the expression of tens of thousands of genes for multiple samples has never been easier. The cost of these experiments has now moved from generating the data to storing and analysing it.

There are many steps involved in analysing an RNA-Seq experiment. Analysing an RNAseq experiment begins with sequencing reads. These are aligned to a reference genome, then the number of reads mapped to each gene can be counted. This results in a table of counts, which is what we perform statistical analyses on in R. While mapping and counting are important and necessary tasks, today we will be starting from the count data and getting stuck into analysis.

First, let’s load all the packages we will need to analyse the data.

# load libraries
library(tidyverse)
library(ttBulk)
library(tidyHeatmap)

GREIN repository of RNA-Seq datasets

In this tutorial, we will learn some R through creating plots to visualise data from an RNA-Seq experiment. RNA-Seq counts file can be obtained from the GREIN platform. GREIN stands for GEO RNA-Seq Experiments Interactive Navigator and provides >6,500 published datasets from GEO that have been uniformly processed. It is available at http://www.ilincs.org/apps/grein/. You can search for a dataset of interest using the GEO code. We obtained the dataset used here using the code GSE60450. GREIN provide QC metrics for the RNA-Seq datasets and both raw and normalised counts. We will use the raw counts here. Generally, the higher the number of counts the more the gene is expressed.

Mouse mammary gland dataset

Here we will perform RNA-Seq analysis using data from a breast cancer research study, from the paper by Fu et al. 2015, GEO code GSE60450. This study examined gene expression in basal and luminal cells from mice at different stages of mammary gland development (virgin, pregnant and lactating). There are 2 samples per group and 6 groups, 12 samples in total.

Reading in the data

Set up an RStudio project specifying the directory where you have saved the /data directory. Open a new script for this workshop File > New File > R Script. Save it as e.g. intro-rnaseq.R.

# read in counts file
counts <- read_csv("data/GSE60450_GeneLevel_Raw_data.csv")

Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_character(),
  gene_symbol = col_character(),
  GSM1480291 = col_double(),
  GSM1480292 = col_double(),
  GSM1480293 = col_double(),
  GSM1480294 = col_double(),
  GSM1480295 = col_double(),
  GSM1480296 = col_double(),
  GSM1480297 = col_double(),
  GSM1480298 = col_double(),
  GSM1480299 = col_double(),
  GSM1480300 = col_double(),
  GSM1480301 = col_double(),
  GSM1480302 = col_double()
)

# read in metadata
sampleinfo <- read_csv("data/GSE60450_filtered_metadata.csv")

Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_character(),
  characteristics = col_character(),
  immunophenotype = col_character(),
  `developmental stage` = col_character()
)

Let’s take a look at the data. You can type the name of the object to view the first few lines and to see how many rows and columns it has.

counts

The counts object contains information about genes (one gene per row), the first column has the Ensembl gene id, the second has the gene symbol and the remaining columns contain information about the number of reads aligning to the gene in each experimental sample. Note the gene counts here are not integers as they’re estimated counts from salmon (see here: https://support.bioconductor.org/p/101156/). There are two replicates for each cell type and time point (detailed sample info can be found in file “GSE60450_series_matrix.txt” from the GEO website). The sampleinfo metadata file contains basic information about the samples that we will need for the analysis today.

First we will convert the counts into long format (tidy format), similar to what we did in the Intro to R session.

# convert to tidy format
counts <- pivot_longer(counts, cols = starts_with("GSM"), names_to = "sample", values_to = "count") 
# take a look
counts

We will next extract just the columns we need, sample, gene_symbol, count. To do this we will use the tidyverse pipe %>%. This ‘pipes’ the output from the command on the left into the command on the right/below. Using the pipe is not essential but it reduces the amount of code we need to write when we have multiple steps (as we’ll see later). It also can make the steps clearer and easier to see. For more details on the pipe see here.

# using pipe
counts <- counts %>% 
  select(sample, gene_symbol, count, X1)
# take a look 
counts

Take a look at the sampleinfo file. The first column “X1” contains the sample ids, the second “characteristics” contains the specific group the sample belongs to (e.g. mammary gland, luminal cells, virgin), the third column “immunophenotype” contains just the cell type (luminal or basal) and the fourth column “developmental stage” contains just the stage (virgin, pregnant or lactating).

sampleinfo

We want to compare the groups in the “characteristics” column however the names are quite long so, similar to what we did in the Intro to R session, we’ll make a column containing shorter group names.

# make column called condition with shorter group names
sampleinfo <- mutate(sampleinfo, condition = case_when(                     str_detect(characteristics, "basal.*virgin") ~  "bvirg",
        str_detect(characteristics, "basal.*preg")  ~  "bpreg",
        str_detect(characteristics, "basal.*lact")  ~  "blact",
        str_detect(characteristics, "luminal.*virgin")  ~  "lvirg",
        str_detect(characteristics, "luminal.*preg")  ~  "lpreg",
        str_detect(characteristics, "luminal.*lact")  ~  "llact"
       ))
sampleinfo

Now we have our counts matrix in the long format we will join it to our sampleinfo so we have information on the samples, what groups they belong to. This is similar to what we did in the Intro to R session.

counts <- full_join(counts, sampleinfo, by = c("sample" = "X1"))
# take a look
counts

Now that we have our data in the format we want we will create a ttBulk object, that we can use to perform differential expression analysis with the ttBulk package. For this we need to specify our counts object and the names of the columns that contain our sample ids, our gene identifiers and our counts. Any other columns in the counts object e.g. our Ensembl gene id “X1” column will remain at the end.

#create a 'tt' object
counts <- ttBulk(counts, sample, gene_symbol, count)
# take a look
counts

Some gene symbols are not unique, they map to more than one gene id. We need to remove this redundancy and we can do that with ttBulk function aggregate_duplicates(). By default it will aggregate duplicate gene symbols summing their counts. Add way to identify what the duplicates are?

# get rid of duplicated gene symbols
counts <- aggregate_duplicates(counts)

Filtering and Normalisation

Before testing for differential expression we need to do some preprocessing of the data. We filter lowly expressed genes and also normalise for differences in sequencing depth and composition between the samples. These steps are explained in more details below.

First we will check how many counts we have for each sample. We can do this quite easily by making a bar plot. This helps us see whether there are any major discrepancies between the samples more easily.

# make barplot of counts
ggplot(data=counts, mapping=aes(x=sample, weight=count, fill=condition)) + 
  geom_bar()

The bar plots show us there are ~20 million counts per sample.

Filtering lowly expressed genes

Genes with very low counts across all libraries provide little evidence for differential expression and they interfere with some of the statistical approximations that are used later in the pipeline. They also add to the multiple testing burden when estimating false discovery rates, reducing power to detect differentially expressed genes. These genes should be filtered out prior to further analysis.

There are a few ways to filter out lowly expressed genes. When there are biological replicates in each group, we favour filtering on a minimum count threshold in the smallest group size. This ensures that a gene will be retained if it is only expressed in one group. In this dataset, we choose to retain genes if they are expressed at a counts-per-million (CPM) above 0.5 in at least two samples (our smallest group size). As a general rule, a good CPM threshold can be chosen by identifying the CPM that corresponds to a count of 10, which in this case is about 0.5 as there are ~20 million counts per sample (20/10 = 0.5 CPM). If the count is any smaller, it is considered to be very low, indicating that the associated gene is not expressed in that sample. Smaller CPM thresholds are usually appropriate for larger libraries. So for this dataset we will filter lowly expressed genes using a cpm_threshold of 0.5 and a prop_threshold (proportion of samples) of 2/12 with ttBulk.

Normalisation for sequencing depth and composition

TMM normalisation is performed to eliminate composition biases between libraries [@robinson2010tmm]. This generates a set of normalisation factors, where the product of these factors and the library sizes defines the effective library size. TMM normalisation (and most scaling normalisation methods) scale relative to one sample.

In the ttBulk package the function scale_abundance() performs the filtering and normalisation to generate normalised counts.

# Normalisation for library size and composition bias (scale counts), 
counts.norm <- counts %>% scale_abundance(
  cpm_threshold = 0.5,
  prop_threshold = 2/12)
# take a look
counts.norm

After we run scale_abundance() we should see some columns have been added to the end of counts. We have a column called filter out low counts that indicates whether the gene has been filtered due to being lowly expressed, FALSE means it wasn’t filtered, TRUE means it was. The count scaled column contains the scaled counts, after the normalisation has been applied.

Maybe add a way to see how many genes we have for each sample, before and after filtering?

We can create density plots to view the distributions of the counts for the samples before and after filtering. This is also a quality check to see if the samples look similar and that none look majorly different. Note we need to take the log2 of the counts and add a small offset (1) to avoid taking log of zero.

# density plot before 
counts.norm %>% 
    ggplot(aes(x=log2(`count scaled` + 1), group=sample, color=condition)) +
    geom_density()

The large peak on the left of the plot shows that a large proportion of genes within each sample are not expressed or lowly-expressed.

Exercise

Adapt the code above to create a density plot of the counts after filtering lowly expressed genes. How does it compare it to the density plot above?

# Solution
counts.norm %>% 
  filter(`filter out low counts` == FALSE) %>%
    ggplot(aes(x=log2(`count scaled`+1), group=sample, color=condition)) +
    geom_density()

We can also create box plots to check the distributions of the counts in the samples. We can add a line through the median to help us see how similar (or not) the distributions are.

# box plot before scaling
counts.norm %>% 
  filter(`filter out low counts` == FALSE) %>%
    ggplot(aes(x=sample, y=log2(count+1), fill=condition)) +
    geom_boxplot() +
  geom_hline(aes(yintercept = median(log2(count+1)), colour = 'red'))

Exercise

Adapt the code above to create box plots that only includes the genes that have not been filtered due to low counts. How does it compare it to the box plots above?

# box plot after filtering
counts.norm %>% 
  filter(`filter out low counts` == FALSE) %>%
    ggplot(aes(x=sample, y=log2(`count scaled`+1), fill=condition)) +
    geom_boxplot() +
  geom_hline(aes(yintercept = median(log2(`count scaled`+1)), colour = 'red'))

Multidimensional scaling plots

By far, one of the most important plots we make when we analyse RNA-Seq data are MDS plots. An MDS plot is a visualisation of a principal components analysis, which determines the greatest sources of variation in the data. A principal components analysis is an example of an unsupervised analysis, where we don’t need to specify the groups. If your experiment is well controlled and has worked well, what we hope to see is that the greatest sources of variation in the data are the treatments/groups we are interested in. It is also an incredibly useful tool for quality control and checking for outliers. We can use the reduce_dimensions() function to calculate the dimensions.

# get MDS dimensions
counts.norm.MDS <-
  counts.norm %>%
  reduce_dimensions(method="MDS", .dims = 2)
# take a look
counts.norm.MDS

Then we can select just the dimensions for the samples.

# get the dimensions with all metadata
MDSdims <- counts.norm.MDS %>%
select(contains("Dim"), sample, immunophenotype, `developmental stage`, condition) %>%
distinct()
# take a look
MDSdims

Next we can plot the MDS dimensions as a scatterplot.

# MDS plot
ggplot(MDSdims, aes(x=Dim1, y=Dim2, colour=sample)) + 
  geom_point()

Exercise

Colour the MDS plot with different metadata variables e.g. immunophenotype. Try using shape= inside the aes(). You can use ?geom_point to check the help page. Discuss what is the greatest source of variation in the data (i.e. what does dimension 1 represent)? What is the second greatest source of variation in the data?

Solution

# MDS plot coloured by cell type
ggplot(MDSdims, aes(x=Dim1, y=Dim2, colour=condition)) + 
  geom_point()

# MDS plot coloured by cell type
ggplot(MDSdims, aes(x=Dim1, y=Dim2, colour=immunophenotype)) + 
  geom_point()

# MDS plot coloured by stage
ggplot(MDSdims, aes(x=Dim1, y=Dim2, colour=`developmental stage`)) + 
  geom_point()

# MDS plot coloured by stage with shape for cell type
ggplot(MDSdims, aes(x=Dim1, y=Dim2, colour=`developmental stage`, shape=immunophenotype)) + 
  geom_point()

Demo more MDS plots (sample swap, batch effects)

Differential expression

Now that we are happy that the data looks good, we can continue to testing for differentially expressed genes. We will use the test_differential_abundance() from ttBulk which uses edgeR to perform the differential expression analysis. We give test_differential_abundance() our ttBulk counts object, our filtering thresholds and a formula (e.g. 0 + condition), specifying the column that contains our groups to be compared. We can also provide the names of the groups we want to compare/contrast to .contrasts (e.g. .contrasts = c(“conditionbpreg - conditionblact”))

# Differential expression with limma-voom (or edgeR)
counts.de <- counts %>%
    test_differential_abundance(
      cpm_threshold = 0.5,
      prop_threshold = 2/12,
      ~ 0 + condition,
      .contrasts = c("conditionbpreg - conditionblact"))

ttBulk says: The design column names are "conditionblact, conditionbpreg, conditionbvirg, conditionllact, conditionlpreg, conditionlvirg" in case you are interested in contrasts
Joining, by = "gene_symbol"
Joining, by = "gene_symbol"

# take a look
counts.de

Now we have columns with our logFC and FDR P values.

Exercise

Perform differential expression for lpreg vs llact

Perform for 2 contrasts, bpreg vs blact and lpreg vs llact at the same time

# maybe include
# remove condition from column names
# counts.de <- counts.de %>% 
#  rename_at(vars(contains("condition")), ~ str_replace(., "_condition", ""))

We can take a look at the top genes by P value. We can use filter() to select only genes with FDR < 0.05, then use distinct() to get the unique values for columns of interest, and arrange() to sort the table by PValue.

counts.de %>%
  filter(FDR < 0.05) %>%
  distinct(gene_symbol, logFC, PValue, FDR) %>%
  arrange(PValue)

We can write out our results to a file that can be loaded into e.g. Excel. write_tsv() will create a tab-separated file.

# save results
write_tsv(counts.de, "my_de_results.tsv")

Plots after testing for DE

Let’s make a few plots to make sure everything looks good and that we haven’t made a mistake in the analysis. Genome-wide plots that are useful for checking are MA plots and volcano plots. We can also use stripcharts and heatmaps to visualise groups of genes.

MA plots

MA plots enable us to visualise amount of expression (logCPM) versus logFC. Highly expressed genes are towards the right of the plot. We can also colour significant genes (e.g. genes with FDR < 0.05)

# maplot, minimal
counts.de %>%
  filter(`filter out low counts` == FALSE) %>%
  ggplot(aes(x=logCPM, y=-logFC, colour=is_de)) +
  geom_point()

# colour by logfc & set up to red, down to blue

Volcano plots

Volcano plots enable us to visualise significance of expression (logCPM) versus logFC. Highly significant genes are towards the top of the plot. We can also colour significant genes (e.g. genes with FDR < 0.05)

# volcanoplot, minimal
counts.de %>%
  filter(`filter out low counts` == FALSE) %>%
  ggplot(aes(x=logFC, y=-log10(PValue), colour=is_de)) +
  geom_point()

# colour by logfc & set up to red, down to blue, highlight some genes

To see how to make more complicated volcano plots, including how to label genes in the plot, see the volcano plot tutorial here. More complicated MA plots could also be made in a similar way.

Stripcharts

In addition to the genome-wide plots already discussed, it is recommended to have a look at the expression levels of the individual samples for the genes of interest, before following up on the DE genes with further lab work. We can use stripcharts and heatmaps to do this. These will help show if expression is consistent amongst replicates in the groups.

For these plots showing expression in the individual samples, we need the normalised counts per gene per sample, so we will first join the scaled counts to the de results.

counts.de <- full_join(counts.de, counts.norm)

Joining, by = c("sample", "gene_symbol", "count", "X1", "characteristics", "immunophenotype", "developmental stage", "condition", "merged transcripts", "filter out low counts")

With stripcharts we can see if replicates tend to group together and how the expression compares to the other groups. Note this is bit more complicated maybe omit. Or break it into steps e.g. get top genes, then plot

# stripcharts
counts.de %>%
    inner_join( (.) %>% distinct(gene_symbol, PValue) %>% arrange(PValue) %>% head(6)) %>%
    ggplot(aes(x = condition, y = log2(`count scaled` + 1), colour = condition)) +
    geom_jitter() +
    facet_wrap(~gene_symbol)

Joining, by = c("gene_symbol", "PValue")

Heatmaps

We can create heatmaps for the most differentially expressed genes. For example we could select the genes with FDR < 0.05 and a logFC change of 4 and make a heatmap of those.

# basic
counts.de %>%
  filter(FDR < 0.05 & abs(logFC) > 4) %>%
  heatmap(
        .horizontal = sample,
        .vertical = gene_symbol,
        .abundance = `count scaled`,
        annotation = c(immunophenotype, `developmental stage`),
        log_transform = TRUE
    )

# customised
counts.de %>%
  filter(FDR < 0.05 & abs(logFC) > 4) %>%
  heatmap(
        .horizontal = sample,
        .vertical = gene_symbol,
        .abundance = `count scaled`,
        annotation = c(immunophenotype, `developmental stage`),
        log_transform = TRUE,
        palette_abundance = c("blue", "white", "red"),
        column_names_gp = gpar(fontsize = 8)
    )

Key Points

• RNA-Seq data can be analysed in a ‘tidy’ way using the packages tidyverse, ttBulk and tidyHeatmap
• Key steps in an RNA-Seq analysis are filtering lowly expressed genes, normalisation for sequencing depth and composition, and testing for differential expression
• MDS plots are very important for examining the quality of the data
• Other useful plots for assessing RNA-Seq data are bar plots, density plots, box plots, MA plots, volcano plots, stripcharts and heatmaps

Further Reading

RNA-Seq analysis is easy as 1-2-3 with limma, Glimma and edgeR RNA-Seq analysis in R